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We simulate the long-range inter-complex electronic energy transfer in Photosystem II - from 
the antenna complex, via a core complex, to the reaction center - using a non-Markovian (ZOFE) 
quantum master equation description that allows us to quantify the electronic coherence involved 
in the energy transfer. We identify the pathways of the energy transfer in the network of coupled 
chromophores, using a description based on excitation probability currents. We investigate how 
the energy transfer depends on the initial excitation - localized, coherent initial excitation versus 
delocalized, incoherent initial excitation - and find that the energy transfer is remarkably robust with 
respect to such strong variations of the initial condition. To explore the importance of vibrationally 
enhanced transfer and to address the question of optimization in the system parameters, we vary 
the strength of the coupling between the electronic and the vibrational degrees of freedom. We 
find that the original parameters lie in a (broad) region that enables optimal transfer efficiency, 
and that the energy transfer appears to be very robust with respect to variations in the vibronic 
coupling. Nevertheless, vibrationally enhanced transfer appears to be crucial to obtain a high 
transfer efficiency. We compare our quantum simulation to a “classical” rate equation based on 
a modified-Redfield/generalized-Forster description that was previously used to simulate energy 
transfer dynamics in the entire Photosystem II complex, and find very good agreement between 
quantum and rate-equation simulation of the overall energy transfer dynamics. 


I. INTRODUCTION 

The initial steps in photosynthesis of green plants occur when Photosystem II (PSII) absorbs light and the nascent 
excitation is used to split water molecules into electrons, hydrogen ions, and molecular oxygen. The ensuing electrons 
and protons drive the subsequent photosynthetic reactions necessary for the production of adenosine tri-phosphate 
(ATP), the “chemical currency” of biological reactions [T]. Greater than 80% of all excitations initiated in PSII will 
result in productive photochemistry, causing the creation of an irreversibly separated electron-hole pair contributing to 
water splitting [5]. The mechanisms of efficient energy transport within PSII have been studied using spectroscopy [3- 
E] and numerical simulations [IMIS], with models based on spectroscopic data and structural information from X-ray 
crystal structures of the complexes [IMS2]- Recent work has focused on energy transfer within PSII supercomplexes, 
which represent the smallest photosynthetic unit that contains all of the relevant proteins for these first steps in 
photosynthesis I1113EI. The largest PSII supercomplex isolated by Gaffarri and co-workers [E] (called G 2 S 2 M 2 ) 
contains 326 pigments bound into an assembly of light harvesting proteins surrounding a PSII reaction center. The 
C 2 S 2 M 2 supercomplex is a dimer, where each of the two monomers contains two antenna complexes (LHCII), three 
minor light harvesting complexes (GP24, GP26, and GP29), two core complexes (CP43, GP47), and one reaction 
center (RC). The arrangement of the complexes is shown in Figure 

Absorption of sunlight by chlorophyll pigments in the antenna complexes of PSII creates electronic excitations that 
are transferred to other pigments located in the same protein complex and from there to pigments in neighboring 
complexes. In this way, the energy is transferred from the antenna complexes to the core complexes and from there to 
the reaction center. Because of these complicated dynamics, involving many degrees of freedom and a relatively large 
number of pigments, numerical simulations for PSII are very challenging. Gonsequently, simulations of energy transfer 
in PSII have been limited to smaller sub-complexes of PSII or to simpler rate-equation descriptions [I5l[l6l[l8| rather 
than the full quantum dynamics simulations that have been performed for smaller light-harvesting systems such as the 
Fenna-Matthews-Olson complex j^MSOj- Recently, in Ref. m, simulations based on a combined modified-Redfield- 
generalized-Forster (MRGF) rate equation approach were carried out for the entire PSII supercomplex, providing 
new insight into the interplay of short-range transfer dynamics inside the individual subcomplexes (proteins) and 
long-range inter-complex dynamics within a rate-equation kinetic analysis. 

In recent years, spectroscopy experiments have provided evidence that there is coherence involved in excitonic energy 
transfer in light-harvesting complexes 1311132 ]. This observation raises questions about the role of such coherence in 
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Figure 1: Structural arrangement of the pigments within the PSII supercomplex, based on the structure determined by Caffarri 
and co-workers [5^. The protein scaffold is removed to show the pigments bound by each protein unit. The pigments that 
belong to the LHCII-CP43-RC subsystem considered in the energy-transport simulations of the present paper are colored: green 
in the LHCIIC monomer, red in CP43, and orange in the reaction center. The other pigments not included in our truncated 
model are gray. 


the energy transfer. In particular, is coherence important for the efficiency of the transfer? Is it relevant for design 
principles that aim to maximize this efficiency? These questions have been addressed in many recent papers |33H45j . 
In a companion paper |46j of the present paper, the role of electronic coherence in electronic excitation energy transfer 
is analyzed in the framework of excitation probability currents. 

In the present paper, we investigate the energy transfer dynamics in PSII by means of a full quantum simulation, 
using a non-Markovian quantum master equation description. This not only allows us to quantify the electronic 
coherence involved in the energy transfer over length scales including several subcomplexes, but also to calculate its 
contribution to the energy transfer in terms of excitation probability currents, using the framework of Ref. |46j . The 
analysis of the probability currents also provides insight into the pathways of the energy transfer in the large network 
of coupled pigments. Using a non-Markovian quantum master equation allows us to account for the non-negligible 
memory times, i.e., the non-Markovianity, associated with both the intra-molecular vibrations of the pigments and the 
vibrations of the protein environment of the pigments to which the electronic degrees of freedom couple. Specifically, 
we investigate the long-range inter-complex energy transfer from the antenna complex via a core complex to the 
reaction center. To this end, we simulate the energy transfer in a subsystem containing a LHCII monomer, the 
core complex CP43, and the reaction center (see Figure]^. This requires incorporation of ~ 30 pigments in the 
simulation. To be able to efficiently treat non-Markovian dynamics for this number of pigments, we use the non- 
Markovian ZOFE quantum master equation that has been developed in recent years m, successfully applied to 
biological light-harvesting complexes, and also tested against exact results EaED]. This approach allows calculations 
for wide parameter ranges of the model. To compare the simulation results of the ZOFE quantum master equation 
and the MRGF rate equation of Ref. m, we analyze here the results of both methods for excitonic energy transfer 
in the subsystem of PSII described above. 

There is an ongoing discussion about whether and how energy transfer in light-harvesting systems depends on 
the initial conditions - in particular, initial excitation in laser spectroscopy experiments can be coherent and local¬ 
ized, whereas in natural initial conditions through absorption of sunlight, delocalized, incoherent states are often 
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assumed [511 [311 SSI 03] ■ Motivated by these discussions, in the present paper we simulate the energy transfer dynam¬ 
ics for very different initial states to ascertain how long-range energy transfer in PSII depends on the initial condition. 
We consider first initial excitation that is localized in the antenna complex, then look at how the energy transfer 
changes when all pigments in the antenna complex, core complex, and the reaction center carry the initial electronic 
excitation, i.e., when the initial state is completely delocalized over all three complexes. We also analyze the role of 
vibrations - in particular, intra-molecular vibrations of the pigments and vibrations of the protein environment - that 
couple to the electronic degrees of freedom and that can enhance the electronic energy transfer through vibrationally 
enhanced transport, i.e., by creating resonances between pigments through vibronic energy levels of the pigments. In 
the present work, we vary the coupling between the electronic and the vibrational degrees of freedom to investigate 
how the energy transfer depends on this coupling to the vibrations. 

The remainder of the paper is structured as follows. In Section |n| we describe the model that we use for the 
simulation of the energy transfer in PSII. The simulated energy transfer dynamics are then shown in Section m 
for an initial state in which the initial excitation is localized in the antenna complex. In Section IV 
to a simulation where the initial excitation is delocalized over all three complexes considered here. 


we compare 
In Section IIVI 


the dependency of the energy transfer on the coupling between the electronic and vibrational degrees of freedom is 
investigated. A summary of the results and concluding remarks are given in Section [^ 


II. MODEL 

We aim to investigate the long-range transport of excitation energy in the Photosystem II supercomplex, from 
the surrounding antenna complexes, through the outer components of the supercomplex, to the reaction center - 
by means of a full quantum simulation. To this end, we model a subsection of the Photosystem II supercomplex 
originally isolated by Croce, et al. [53] and subsequently used as a foundation for a structure based excitation energy 
transport model to explain the measured fluorescence lifetimes [5]. Figure shows the largest PSII supercomplex 
previously modeled mm- The colored pigments in LHCII, CP43, and the reaction center shown in Figure [^represent 
the subsystem that we model in this work: it contains 33 chromophores distributed between a LHCII monomer (14 
pigments), CP43 (13 pigments), and a truncated reation center (6 pigments). 


A. System parameters 

1. Electronic degrees of freedom 

Our model contains 33 chromophores in total, each of which we model as having two electronic states separated by 
the local excitation energy, referred to here as the ’site energy’, which varies across different pigments depending upon 
their interactions with the local protein environment. We further describe the interaction between the chromophores 
by electronic matrix elements that couple electronic excitation of a chromophore to electronic excitation of other 
chromophores. In the subspace of single excitations that is relevant for energy transport, this is described by the 
electronic Hamiltonian, 


Helec = y^£n| n)(n I + Vn,Tn\n){m\ (1) 

n nm 

containing N site energies (Sn) and the N x N coupling matrix V. Here, |n) is the state in which only pigment 
n is excited and all others are in the ground state. The construction of an electronic Hamiltonian for the PSII 
supercomplex has been described previously m- The couplings between pigments contained within the same protein 
have been extracted from the literature for each complex |16] , where they were constructed to reproduce the available 
linear and non-linear spectroscopy data for ensembles of the isolated pigment-protein complex. The coupling between 
pigments in different proteins has been calculated assuming a dipole-dipole coupling m- All energies and coupling 
parameters are taken from Ref. [16], and are shown in Figure]^ together with the structure of the LHCH-CP43-RC 
truncated supercomplex. 


2. Coupling to vibrations 


The electronic states of each chromophore are coupled to a large collection of vibrational modes that represent 
intra-molecular vibrational modes as well as vibrational modes of the surrounding protein scaffold. In the extreme 
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Figure 2: Electronic part of the model used for the simulation of energy transport from the LHCII antenna complex via the 
CP43 core complex to the reaction center (only monomer C of LHCII is included). This subsystem of PSII contains in total 
33 pigments. Also included are the two radical pair states RPl and RP2 of Ref. m with the characteristic transfer times for 
the description of the charge transfer. Figure A: spatial and energetic structure. The spheres are located at the positions of 
the pigments and their sizes show the relative electronic transition energies of the pigments. The numbering of the pigments 
starts with pigment 0 in LHCII that has the highest energy (15,890 cm“^), then the number increases with decreasing energy 
within each protein; pigment 32 has the lowest energy (14,740 cm“^). The connection lines between the spheres represent 
the electronic interactions between the pigments; their thickness is proportional to the strength of the respective interaction 
(absolute values); the maximum interaction strength is 150 cm“^. Transfer of population between pigments 27, 28, and 32 and 
RPl, and from RPl to RP2 is depicted by arrows. In addition to the numbering of the pigments used in the present work, 
the common names of the pigments used in Ref. [16] and references therein are given in the listing. Figure B: corresponding 
Hamiltonian matrix in the site (pigment) basis, where the relative transition energies of Figure A are on the diagonal and the 
absolute values of the interactions are the off-diagonal elements. For better visibility, the interactions are increased by a 
factor of ten with respect to their true values. The coloring at the edges indicates to which protein the sites belong: 
Green: LHCII C. Red: CP43. Orange: RC. Black: RPl state. Blue: RP2 state. Gray: ground state. All parameters for this 
model are taken from Ref. |16| . 


limit of low-frequency vibrations, the long-time-scale protein conformation dynamics give rise to an inhomogeneous 
distribution of energies for each chromophore within any ensemble measurement. In contrast to Ref. m, however, 
where ensemble averaging over static disorder in the site energies is performed, in the present work we do not 
describe ensemble averaging and instead use the electronic Hamiltonian with average site energies for our simulations. 
While ensemble averaging can be critical for comparing to experimental measurements, in this paper we investigate 
quantities that are primarily focused on the fundamental process of excitation energy transport within an individual 
PSII supercomplex. Here we may take the average energies instead of random energies of a single realization to create 
a “best-case scenario” for transport efficiency and emergence of coherence. 

We describe here the intra-molecular vibrational modes of the pigments together with the vibrational modes of the 
protein environment by the single vibrational Hamiltonian 

Hvih = EE (2) 

n X 

where we approximate the vibrational modes to be harmonic. Here, Unx is the annihilation operator of mode A with 
frequency LOnX that belongs to a pigment n. (Here and throughout the paper we set h = 1). In this description, 
each pigment has its own separate bath of vibrational modes that does not directly couple to the modes of another 
pigment. Indirectly, however, such a coupling can come about through the electronic dipole-dipole interaction between 
two pigments. 

We assume that electronic excitation of a pigment couples linearly to its own vibrational modes, such that the 
overall coupling between the electronic degrees of freedom and the vibrations is given by 

.I^elec—vib — ^ ^ ^ ^ ^nxi^liX 4” ^nx): 

X 


n 


( 3 ) 
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where the system operators are the projectors Ln = — |n)(n|. The coupling constants Kn\ that describe the 
strength of the coupling of electronic excitation of pigment n to vibrational mode A of this pigment, can be expressed 
through the vibrational spectral density of pigment n 


^ ^ I^uaI S(cU COfix)- 


Assuming that the spectral density is a continuous function, this leads to the vibrational correlation function 


«n(r) = / dw Jr 

Jo 


(w) coth ( vTTf;) cosIwt) — 1 sinlwr) 
L \2kl / 


which depends on the temperature T. The Fourier transform 


/ OO 

dr e"‘^'"an{T), 

-OO 


( 4 ) 

( 5 ) 

( 6 ) 


which we shall refer to as the ’vibrational correlation spectrum’, will be useful later on for application of the ZOFE 
master equation. It can also be directly calculated from the spectral density Jniui) (or vice versa) via the relation 


C„(a;) = 2 


1 + 




with the anti-symmetrized spectral density 


Jn{^) 


Jn{oj) ; w > 0 

'Ai( 5 cu 0 


( 7 ) 

( 8 ) 


Peaks in C„{uj) that are narrow compared to the relevant energies of the electronic degrees of freedom will lead 
to non-Markovian dynamics. When we apply the ZOFE master equation below, we will fit Cn{oj) using a sum of 
Lorentzians. 


3. Charge transfer in the reaction center 

Up to this point, we have described a system where electronic excitations of individual pigments are coupled both 
to other pigments and to vibrational modes. Two other key features need to be incorporated into our model to make 
it a reasonable model of the initial stages of photosynthesis in PSII. These are, i) energy capture via charge separation 
at the reaction center, and ii) loss of excitation via non-radiative or fluorescence processes. 

We include these dissipative processes by extending our description of the electronic degrees of freedom to include 
the phenomenological radical pair states (|RP1) and |RP2)), as described in Ref. [IS], and the ground state |0), 
where all pigments are in their ground electronic states. 

In this truncated supercomplex model, energy is transferred from the excited states of three of the pigments in 
the reaction center to the first radical pair state RPI. On a slower time scale, some of this excitation is transported 
back to the excited states of the three pigments in the reaction center. From RPI, the energy is irreversibly and 
more slowly transferred to the second radical pair state RP2. The latter is the last step in our model and describes 
the irreversible trapping. The characteristic transfer times that are found in Ref. m are 0.64 ps from the excited 
states of the RC pigments to RPI, 160 ps from RPI back to the excited states of the RC pigments, and 520 ps for 
the irreversible transfer from RPI to RP2. Electronic excitations present on any pigment are assumed to have equal 
probability of being lost through either fluorescent or non-radiative decay mechanisms, with a characteristic decay 
time of 1.78 ns (combining the decay times of 2 ns for non-radiative decay and 16 ns for fluorescence used in Ref. [lOjb 

The Hamiltonian term describing the radical pair states and the ground state is given by 

77rp+gs = egsI 0) (0 I + eRPi I RPI) (RPI I + eRP21 RP2) (RP2 |, (9) 

where egs is the ground state energy, and £rpi and £rp 2 are the energies of the RPI and RP2 state. 

We shall describe in detail the dynamical model employed to treat excitation transfer into either the radical pair 
or ground state of the pigment in Section |IIB 3[ 
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B. Excitation energy transport 


Based on the previous section (Section II A), the total Hamiltonian of our model is then 

Htot ~ ^e\ec A -^RP+GS -^elec—vib d” -^vib- 


( 10 ) 


Since this Hamiltonian comprises a large number of degrees of freedom, direct calculation of the excitation energy 
transport by explicitly solving the full quantum dynamics is not realistic. Therefore, we use instead the framework 
of the description of open quantum systems, and divide the total Hamiltonian into three components: a system 
(Hsys), an environment (i?env)j and the interaction between system and environment (i?sys-env)- The system part 
Hsys should contain the quantities relevant to the energy transport, that is, most importantly the probabilities of 
electronic excitation of the different pigments and the excitation of the radical pair states in the reaction center. On 
the other hand, we want to keep Hgys as small as possible to keep the calculation numerically manageable. Therefore, 
we choose this to consist only of the electronic degrees of freedom of the pigments and the radical pair states. The 
vibrations are then treated as the environment. That is, we have 


^sys — ^elec d” I^RP+GSi ^env — ^vihj and i?sys — env — ^elec—vib- (11) 

The excitation energy transport and the electronic coherence can be extracted from the reduced density matrix of 
the system part, i.e., 

p{t) = Trenv Ptot{t), (12) 

which is obtained from the total density matrix by tracing out the degrees of freedom associated with the environment, 
i.e., the vibrations. In the basis of the states | n) of the local excitations of the pigments, the diagonal of the reduced 
density matrix gives the populations of the excited electronic states of the pigments, the radical pair states, and the 
ground state. The off-diagonal elements are the electronic coherences between the pigments, which constitute a major 
focus of the present work. 


1. ZOFE master equation 


For the simulations of the dynamics of energy transport and coherence we use the ZOFE master equation, which 
allows us to take into account the non-Markovian effects in the coupling between the electronic degrees of freedom 
and the vibrational environment. The ZOFE master equation is given by HSIST] 

dtp{i^ = -I [ilsys, p\+Yl {LnpO^^^^ + O^^'^pLi - LiO^^^p - (13) 

n 

where the auxiliary operator 

[ dsan{t - (14) 

captures the non-Markovian effects, and with initial condition 

0i^\t,t) = L„. (15) 

Depending on the form of the environment correlation function an{t — s), a closed evolution equation for the auxiliary 
operator (f) can be obtained [13 |13 . For an environment correlation function that is a sum of damped 

oscillating terms, 

a„(t - s) = (16) 

3 


such a closed evolution equation can be found [29]. Writing the operator Oo"^(f) as the sum the 

closed evolution equation to obtain the operators {t) is given by [Hj 


dtO^Q^\t) = TnjLn - {iflnj +7nj)Oo”'^^ 


iff 




(17) 


with initial condition = 0) = 0 and = L„. In the present work, we use these evolution equations 

with coupling operators = — | n)(n | for each pigment n. 
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Figure 3: Spectral densities of intramolecular and protein vibrations that the different pigments in the different proteins couple 
to. The thin gray bars at the upper edge show where all transition energies of the electronic system Hamiltonian are, that is, 
all energy differences between the eigenenergies of the Hamiltonian. Figure A: The blue and green narrow peaks are the two 
different spectral densities of the Chi B (blue) and Chi A (green) pigments in the LHCH antenna, respectively, from Ref. |16) 
and references therein. The broad blue and green peaks are the corresponding fitted spectral densities that are used in the 
ZOFE simulation. The green bars at the upper edge show the transition energies between the exciton states localized in the 
LHCH antenna. Figure B: The thin red and orange curves show the spectral densities of the pigments in CP43 (red) and RC 
(orange), respectively, from Ref. |16| and references therein. The thick curves are the corresponding fitted spectral densities for 
the ZOFE simulation. Red, magenta, and orange bars at the upper edge indicate the transition energies between exciton states 
localized in CP43, localized in both, CP43 and RC, and localized in RC, respectively. We note that the spectral densities in 
Figure A are much larger than the ones in Figure B (note range of y-axis). 


2. Approximate environment spectral density for ZOFE master equation 


The form of a„(t — s) in Equation (16) corresponds to a sum of Lorentzians centered at frequencies Flnj and with 


weights Tnj and widths in the environment correlation spectrum C„(w) (the Fourier transform of Q!„(t — s)): 


C„(a;) = -^r„, 

TT • * 


Inj 


(W - 0„j )2 + qT ' 


(18) 


We use this expression to fit the environment correlation spectra Cn{uj) of the pigments to experiment-based spectra, 
in order to obtain the parameters Flnj, ^nj^ and ^nj that we need when we solve the evolution equation for the 
auxiliary operator Eq. The resulting spectral densities Jniyj) (Eq. @) that we use in the ZOFE simulation 

are shown in Figure]^ together with the original, measurement-based spectral densities of Ref. [H]. As one can see 
in Figure [^, the htted spectral density that we use for the LHCH pigments in the ZOFE simulation is very broad 
compared to the narrow peaks of the original spectral density. We approximate the original spectral density with 
these broad peaks, because the ZOFE approximation that we use in the simulation is quite limited in the degree of 
non-Markovianity that it can handle. Thus, the memory time of the vibrations can not be too long compared to the 
time scale of the electronic dynamics. Narrow peaks would result in a long memory time, so we are limited to broad 
peaks that capture only the general features of the coupling to the vibrations. 

The approximate spectral densites in Figure also neglect the high-energy contributions of the original spectral 
density around and above 1,500 cm“^. Since these peaks are off-resonant with the system energies, which are marked 
by the bars at the upper edge, we assume that they do not have a large influence on the system dynamics and it is 
reasonable to neglect them. In Ref. [55], it was shown that the influence of off-resonant parts of the spectral density 
on excitation energy transport is indeed negligible. 

In Figure]^, the spectral densities for CP43 and the RC that we use for the ZOFE simulation are shown. They 
are only in rough agreement with the original, measurement-based spectral densities from Ref. m- This difference 
arises because these spectral densities are obtained not by directly fitting the spectral densities Jn(w), but by fitting 
the environment correlation spectra Cn{i-o) with Lorentzians, according to Equation (18). To summarize, we apply 
the following procedure: 


1. Calculate the environment correlation spectra C„(a;) from the experimentally determined spectral densities 
Jn{oj) from Ref. [T5], using Equation ([^. For this step, we need to assume a temperature that enters in C'„(a;). 
Since we want to simulate energy transport at room temperature, we choose T = 300 K. 

















2. Fit these Cn(uj) with Lorentzians according to Equation (18 1 , to determine the parameters for the evolution 
equations ( [l7| for the ZOFE auxiliary operators. In this work, we use only two Lorentzians for the LHCII 
pigments and one Lorentzian for the CP43 and RC pigments, to keep the corresponding number of evolution 
equatons small (see index j in Eq. 0)- 


As a consistency check, we convert these fitted sums of Lorentzians Cn(uj) back to spectral densities Jn{uj) to compare 
with the original spectral densities, as is done in Figure]^ 


3. ZOFE master equation: excitation loss and radical pair states 


To simulate the transfer of energy to the radical pair states in the reaction center, and to take radiative and non- 
radiative decay of the excited electronic states of the pigments into account, we add Lindblad-master-equation terms 
to the ZOFE master equation. An analogous treatment of radiative decay and trapping in the reaction center is 
employed in Ref. |28j , where a non-Markovian master equation is also extended by corresponding Lindblad terms, for 
the simulation of energy transport in the Fenna-Matthews-Olson complex. (See also Refs. [36l|52] for related Lindblad 
models for trapping and decay.) 

This Markovian Lindblad description of the trapping in the reaction center is a rough approximation of the full 
charge separation dynamics, which we use here in order to create a minimal model to capture the primary dynamics. 
We assume that the characteristic transfer times to the radical pair states that are found in Ref. [16] give a reason¬ 
able approximation and therefore describe the trapping part of our simulation with Lindblad terms based on these 
characteristic transfer times. 

In the Markov limit, where the environment correlation function a(r) decays fast enough compared to all relevant 
time scales of the dynamics, the ZOFE master equation itself becomes a Lindblad master equation |47j. It is therefore 
consistent to add Lindblad terms to the ZOFE master equation, since they are equivalent to additional ZOFE master 
equation terms with fast decaying correlation functions. In the Markov limit ai^r) = 7 ^ 6{ t), th e a uxiliary operator 
of the ZOFE master equation becomes the constant operator OQ*^(t) = ^jiLi (see Eqs. (141 and (15)) [U]- (Here, we 
have replaced the index n, which is used above to refer to the pigments, by a more general index i, to refer to the 
respective process described by a system operator Li and a corresponding coupling constant 7 ^; e.g. relaxation from 
a state | / ) to a state |t) is described by L/t = \t){f\.) Inserting this constant auxiliary operator into the ZOFE 
master equation (13) gives the well-known Lindblad master equation 


dtp{t) = -i [H^ys, P] + X! ( 


L^pLl - ^LlLip- ^pLlLi 


1 




1 


rV: 


(19) 


with coupling (Lindblad) operators Li and corresponding coupling constants 7 ^. We use such a Lindblad description 
for the following processes: 

1. Radiative and non-radiative decay of the electronic excitation of all pigments (n = 0 ... 32) is described through 

Lindblad operators = | 0) (n | and the corresponding coupling constant is given by the decay 

rate. Here, | 0) denotes the electronic ground state. The decay rate is chosen such that it corresponds to a 
characteristic decay time of 1.78 ns (combining the two decay times of 2 ns for non-radiative decay and 16 ns 
for fluorescent decay assumed in Ref. m)- 

2 . Transfer of excitation from the reaction center pigments 27, 28, and 32 to radical pair state RPl through 
L^^^ = |RPl)(n| with rates 7 ^^^ (where n = 27,28,32). From all three pigments the transfer to RPl is 
assumed to occur with the same rate, which is set at the transfer time of 0.64 ps found in Ref. [16]. 

3. Transfer from RPl back to the reaction center pigments 27, 28, 32 through L^^ = \ n)(RPl | with rates 7 ^^. 
It is assumed that equal amounts of population are transferred back to the three pigments, so that each of the 
three rates is one third of the overall rate corresponding to the characteristic transfer time of 160 ps found in 

Ref. [IB] . 

4. Irreversible transfer from RPl to radical pair state RP2 through L^^^ = \ RP2)(RP11 with a rate 7 ^^^ xhis 
rate is chosen according to the characteristic transfer time of 520 ps found in Ref. m- 
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For each of these processes, a corresponding Lindblad term is added to the ZOFE master equation, so that the 
complete master equation we solve in our simulations is 


dtp{.t) 


-i p] + ^ - pOlP 

n—0 

n V / 


( 20 ) 


where the Lindblad operators Li and corresponding coupling constants (rates) 7 ^ in the Lindblad term in the second 
line of the equation belong to the respective processes that are listed above, and the index i of the sum refers to the 
individual processes. In the ZOFE term in the first line, the ZOFE coupling operators L„ = — |n)(n| couple only 
electronic excitation of the pigments to the non-Markovian vibrational environment. The radical pair states | RPl) 
and I RP2) are not coupled to the non-Markovian vibrations. 

The system density matrix pit) now has components for the states | n), in which pigment n is excited, as well as 
components for the radical pair states |RP1) and | RP2) and for the ground state | 0). pit) is obtained by solving 
this master equation together with the evolution equation for the ZOFE auxiliary operators. For the truncated 
supercomplex considered here, pit) is then a 36 dimensional square matrix. The system Hamiltonian does not contain 
coupling to or between RPl, RP2, and the ground state; all coupling to and between these states arises through the 
Lindblad terms. 


4- Population currents and the contribution of coherence 

The transfer of excitation energy between the pigments can be expressed in terms of currents of excitation probability 
between the pigments. As shown in a companion paper |46| . this population current description straightforwardly 
identifies the contribution of coherence to the energy transfer, and this description also clearly separates the respective 
contributions of unitary dynamics, dephasing, and relaxation to the population currents between the pigments |46) . 
Another advantage of the population-current description is that it clearly shows the individual pathways that the 
energy transfer takes between the pigments. By integrating the currents over time, it is then possible to see the 
respective amounts of population that have been transported via each pathway. This insight is useful to identify and 
analyze the functionality of the light-harvesting supercomplex and its constituent complexes and pigments. 

Because the total electronic excitation probability is conserved across the pigments - that is, the reduced electronic 
density matrix fulfills J^nPnnit) = 1 at all times - a continuity equation [15] 

dtPnnit) = ^ jmn(t) (21) 

m^n 


holds, where jmn{t) is the (net) population current at time t that transports population (i.e., excitation probability) 
from a pigment m to another pigment n. When jmnit) is positive, excitation is transported from pigment m to 
pigment n; when it is negative, the energy transfer goes in the other direction, i.e., from pigment n to pigment m. 
Based on the continuity equation (21), for quantum master equations of the form used in the present work Eq. (20), 


it is shown in Ref. [46j that the individual population currents jmn it) between the pigments can be calculated through 


jr.nit) = .7'r‘^’'"(t) + jiT^^it) + f^':-it) 


mn it)) + 0 


(TmnPmm(^) ^nmPnn it)): 


( 22 ) 


where and are the respective contributions of unitary dynamics, dephasing, and electronic 

relaxation, which are specified in the second line of the equation. Here, Pmnit) are the coherences between the sites, 
Pnnit) are the populations of the sites, and Vmn are the inter-site couplings and is the rate of electronic relaxation 
from a site (or state) m to a site n. 

When we have the density matrix pit) at a time t from a numerical simulation as described above, we can use 
Equation (22) to calculate the currents between the sites. As described in detail in Ref. [46], the first term in 
Equation (22) stems from the unitary part of the dynamics - that is, in our model of PSH from the unitary transfer 


driven by the electronic coupling between the pigments. The second term stems from dephasing through the coupling 
to the vibrations and gives zero, that is, the dephasing does not contribute to the population currents. The third term 
stems from electronic relaxation processes, that is, it describes the transfer of population to the radical pair states 
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enabling the charge separation in the reaction center, as well as the radiative and non-radiative decay to the ground 
state. 


We see in Equation (22) that the unitary part of the currents comes entirely from the coherence between the sites 


- more precisely the imaginary part of the coherences. The relaxation part of the currents that describe the transfer 
to the RP states, on the other hand, comes entirely from the populations of the sites, as can be seen in the third term 
in Equation (221. This means that the energy transfer from LHCII to CP43 and to the RC derives entirely from the 


coherence between the sites, whereas the transfer to the RP states is entirely through the population. 

When the coupling to the vibrations is strong compared to the electronic inter-pigment coupling, the resulting 
strong dephasing destroys a large part of the coherence between the pigments. We see from Equation (22) that this 


will decrease the energy transfer between the pigments, because the unitary contribution to the population current, 
which in our model amounts to all the energy transfer between the pigments, decreases since it is entirely described 
by the coherence. This suppression of the transport at strong coupling to the vibrations is well known and often 
explained in the picture of the quantum Zeno effect |53j : the electronic excitation interacts with the vibrations - 
is “measured” by the vibrations - so strongly (i.e., at such a high rate) that the excitation is forced to stay in the 
original state and cannot move anymore. 

To obtain the overall net currents between the different proteins, we first calculate the currents jmn(t) between 
individual sites m and n of the different proteins and then sum these currents up. That is, the current JAsit) between 
subcomplex A and a subcomplex B is given by 




(23) 


m^A nGB 


where when the current JAB{t) is positive, there is a net flow from A to B, and when JAB(t) is negative, there is a 
net flow from B to A. 


III. ENERGY TRANSFER: SIMULATION RESULTS 

In order to allow for a full quantum simulation of energy transfer, we construct a minimal model of light harvesting in 
PSII. In their natural environment (the thylakoid membrane), PSII supercomplexes are surrounded by light harvesting 
proteins (LHCII). Excitation energy transfer from the surrounding LHCII pigments into the supercomplex occurs via 
the light harvesting protein mechanically bound to the core complex (LHCH-S and LHCH-M). Our minimal model of 
the process of excitation energy transfer from the surrounding light harvesting proteins to the reactions center contains 
one LHCII monomer, the complete CP43 protein and a slightly truncated RC. We use the ZOFE master equation 
to calculate the time-dependent density matrix in the space of the electronic states of the pigments, as described in 
the previous section. This enables us to analyze the transfer of the electronic excitation energy and the electronic 
coherence between the pigments. 

It is important to differentiate between the characteristics of excitation energy transfer and the functional behavior 
of the light harvesting apparatus. In low light conditions, plant growth depends on efficient light harvesting - creating 
as many productive photochemical events in the RC as possible per photon absorbed [S3]. As such, in this paper, 
we will equate the fraction of excitations that cause productive photochemistry (i.e. reach the RP2 state) with the 
function of the light harvesting antenna. A decrease in the fraction of excitations reaching RP2, then, is a decline in 
the functional behavior of the photosynthetic assembly. It follows that changes in the precise nature of the excitation 
dynamics do not necessarily result in changes to the functional behavior. Two different excitation dynamics can 
still give rise to equivalent quantities of excitation causing irreversible electron loss (reaching the RP2 state). As 
a result, to show that certain terms in an equation of motion are important to photosynthesis requires both that 
they change the excitation dynamics and that they increase (or decrease) the fraction of excitation that reaches 
RP2, the latter being a measure of how much they influence the function of photosynthesis. In particular, we use 
the amount of population that reaches the second radical pair state RP2 in the reaction center within 1 ns as a 
measure for efhciency of the transport, since this state undergoes irreversible charge separation in our model. (We 
note, however, that fluorescence-lifetime measurements on PSII are only sensitive to the total population remaining 
in chlorophyll excitation. As a result, the experiment cannot directly differentiate the kinetics of electron transport, 
and the assignments to RPl and RP2 must be considered phenomenological [16].) 

In our first study, we describe excitation energy transfer when initial excitation occurs in a pigment-protein complex 
outside the supercomplex and subsequent energy transfer occurs into the core complex via the attached LHCII trimer. 
We therefore initiate excitation on two chlorophyll A molecules in the LHCII monomer (sites 7 and 10) to represent 
excitation that is entering the system from the neighboring LHCII monomers, since these sites are assumed to have a 
large rate of transfer with the neighboring LHCII monomers [T5]. Each of the sites 7 and 10 carry half of the initial 
population, and we choose the initial state to not contain any coherence between the sites. Thus, all off-diagonal 
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Figure 4: Energy transfer dynamics calculated with the ZOFE master equation, for the electronic model and vibrational 
environment described in Section [IT] Initially, only sites 7 and 10 are excited, with the same amount of population on both 
pigments, and without initial coherence in the site basis. Figure A: populations in the different proteins LHCllC, CP43, and 
in the RC, and of the two radical pair (RP) states RPl and RP2. For comparison, the results from the generalized-Forster- 
modified-Redheld calculation of Ref. [1^ (dashed lines) and of a pure Forster calculation (dotted lines) are shown. (Up to 0.1 
ps, the time axis is linear; after that it is logarithmic). Figure B: po pulati on currents between the proteins and RP states, 
calculated from the time-dependent density matrix using Equations \22\ and (23l. The coherence and population contributions 
to the currents are shown. (The currents have the directions specified in the legend when they are positive, and the reverse 
direction when they are negative). Figure C: total amount of coherence C(p) = m^n |Pn,m|. Solid line: in site basis. 
Dashed: exciton basis. Dotted: in the domain-exciton basis used for the Forster-Redfield simulation of Ref. [16] . Figure D: 
elements of the density matrix in the site basis integrated over time, up to 1 ns. On the upper triangle, the absolute values of 
the imaginary parts of the coherences are shown, that is, J |Imag(p„m)|df. On the lower triangle, the real parts J \R,e{pnm)\dt 
are shown, and on the diagonal the populations J pnndt are shown. (The RPl and RP2 populations are scaled down to 1% of 
their true values. As in Figurej^ the coloring at the edges indicates to which protein the sites belong, and the black lines mark 
the corresponding blocks of the matrix). 


elements of the initial density matrix are zero in the site basis. The dynamics following delocalized initial excitation 
will be discussed in Section ITV Al 


A. Population 

To investigate the (long-range) excitation energy transfer between the complexes and the transfer of energy to the 
charge separated states (RPl and RP2), we consider the population of the electronic excitation in the individual 
complexes and the population of the RP states over time. Figure shows these population dynamics. The solid 
lines in Figure]^ are the populations on each complex calculated with the ZOFE master equation. The population 
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of site 7 in LHCII is shown separately (thin green line). Nearly all excitation (population) is transferred from LHCII 
through CP43 and RC to the first radical pair state RPl within about 200 ps. Irreversible transfer from RPl to RP2 
takes much longer because of the much weaker coupling between RPl and RP2. After 1 ns, 80% of the population 
has been transferred to RP2. Within this time, only a small fraction of less than 5% is lost through radiative and 
non-radiative decay, in keeping with the relatively slow dynamics of excitation loss in a photosynthetic apparatus. As 
we shall see in the following discussion, this separation of timescales between excitation energy reaching the RC and 
the relatively slow loss mechanisms of fluorescence and non-radiative decay is a key feature of efficient light harvesting 
in PSII. 

The dashed lines in Figure show the population calculated with the modified-Redfield/generalized-Forster 
(MRGF) method of Ref. [16]. In the MRGF calculations, pigments are grouped together into “domains”. The 
assignment of pigments to different domains is determined by a threshold electronic coupling and the distribution of 
site energies. Within these domains, the dynamics are treated with modified Redheld theory. Between the domains, 
the transport is calculated based on transfer rates between the exciton states of one domain (electronic eigenstates of 
the domain) and the exciton states of another domain. These rates are determined using generalized Forster theory, 
i.e., employing the overlaps between emission and absorption spectra of the different domains m- 

As we can see in Figure]^, the population dynamics from the MRGF calculation - especially the RP2 population 
- agree fairly well with that of the ZOFE simulation, even though the former is a much simpler classical rate equation 
calculation. This is consistent with the findings of Ref. |36| that classical rate equations are adequate for the calculation 
of transport efficiencies when they are properly derived from the full quantum description. 

In Figure]^, we also show the result from a pure Forster calculation (dotted lines) for comparison, since this simple, 
perturbative method is often used to calculate energy transfer in biological systems. The pure Forster calculation 
is based on rates obtained from the overlaps between emission and absorption spectra of each pair of pigments (i.e. 
donor and acceptor pigments) |551156j . This simple approximate method can give reasonable results in a regime where 
the coupling between donor and acceptor pigments is weak compared to the coupling to the vibrations. However, 
since this condition is not met for many pigment pairs in PSII, especially inside the more strongly coupled domains, 
we expect the pure Forster calculation to be less accurate than the ZOFE calculation as well as the MRGF calculation 
of Ref. [ini . 

In Figure]^ we do indeed see a signihcant deviation of the pure Forster calculation (dotted lines) from the other 
simulations. Furthermore, the transport efficiency (population reaching RP2 within 1 ns) calculated with pure Forster 
is lower. 

The observation that the MRGF and ZOFE calculations are in close agreement with each other, while the pure 
Eorster calculation gives a lower transport efficiency, supports the key role of exciton delocalization in excitation 
transport through PSII. The dynamic effect is sometimes referred to as supertransfer, since it occurs when strongly 
coupled pigments allow an excitation to coherently delocalize within the donor domain, thereby enhancing the rate 
of transport by up to a factor equal to the number of pigments in the donor domain (see e.g. Ref. [39j and references 
therein). This enhancing effect of supertransfer due to excitonic delocalization within a subcomplex or domain 
is naturally included in the MRGF and ZOFE calculations, but it is not taken into account in the pure Forster 
description, which thus leads to a lower transport efficiency. 


B. Population currents 


Figure]^ shows the population currents j(t) between the proteins and from the reaction center to the radical pair 
states. These currents are calculated from the time-dependent density matrix using Equations (221 and (23). 

From Equation (221, we know that the energy transfer from LHGII to CP43 and on to the RG is entirely through 
the coherence between the pigments, and the transfer to the RP states is entirely through the population. This 
is combined by Figure j^, where the solid lines are the total currents, and the dashed and dotted lines show the 
coherence and population contribution, respectively. 

We see that all currents between complexes are positive over the entire time range. That means that the net flow 
of energy is directed towards the RP2 state at all times and no temporary intercomplex backflow occurs. Because 
of the high rate of transport between the RG and the RPl state, in Figure]^ the current from RG to RPl almost 
completely follows the current from CP43 to the RG; there is only a very small delay of the current to RPl. That 
means the population is almost passed right through the RG to RPl. The current to the RP2 state has a smaller 
magnitude but stretches over a longer time range. Therefore the overall population that is transferred to RP2 during 
this time span, which is simply the current integrated over time, is not much smaller than the population transferred 
to the RG, namely 80% versus > 95%. As expected, the direct current between LHCII and the RC (cyan curve in 
Figure]^) is zero, due to the large distance between the LHCII and RC pigments. 

To see the individual pathways of the energy transfer, we integrate the population currents between the individual 
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Figure 5: Population currents between the individual chromophores, integrated over time, that is, J jmn(t) dt. Integrated up to 
a time of 1 ns, and for the same simulation as Figure]^ As in Figure]^ initially only pigments 7 and 10 are excited. The arrows 
show the directions and their thickness the relative magnitude of the integrated currents. Only integrated currents larger 0.03 
are shown. The gray arrows show the currents inside the proteins and the cyan arrows the ones between the proteins. The 
relative positions of the pigments in 3D are obtained from Ref. [T5]. (The labeling of the pigments is the same as in Fig.|^) 


pigments over time. The result is shown in Figure We see in Figure that there is a complex network of pathways 
rather than just one or two dominant pathways. Previous simulations using the MRGF model have shown that 
transport through PSII does not occur by a single, obligate pathway |16j . Our current simulations demonstrate that 
there are multiple transport pathways from LHCII to CP43, and from CP43 to the RC. However, particularly within 
complexes, there are a few pathways that transport a particularly large amount of excitation energy, identified in 
Figure 

Remarkably, in Figure we see that there are a number of relatively strong currents that run in a circular fashion, 
where energy is transported in a loop between pigments. For instance, there is a loop in LHCII from pigment 7 to 13 
to 2 to 3, and back to 7. In CP43, there are several such loops. For instance, from 25 to 26 to 23, and back to 25. 
And there is a longer one from 25 to 26 to 18 to 22 to 23, and back to 25. There is also one from 16 to 21 to 24, and 
back to 16. 

In the reaction center, there is a strong circular current from 32 to 27 to 28, and back to 32. We also see in Figurej^ 
that in the reaction center overall the dominant net flow to the RPl state is via pigment 30 and 32. Even though 
pigments 27 and 28 are part of the strong circular current, they do not appear to contribute significantly to the net 
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transport to RPl. There is even a net backflow from the RPl state to the pigments 27 and 28, which acts to decrease 
the efficiency of the transport. Pigments 27 and 28 receive a small amount of population from CP43 pigments, but the 
backflow from RPl to these two pigments is larger and therefore seems to outweigh their contribution to the transport 
to RPl. In the light of these pathways, it appears that the efficient excitation energy transfer to and in the RC relies 
mainly on the two pigments 30 and 32, and that therefore other reaction center pigments may rather be needed for 
their role in the electron transport (described by the RP states in our model) than for the excitation transfer. To 
test this hypothesis, we ran additional simulations where - except for the pigments 30 and 32 - we decoupled all RC 
pigments from the rest of the pigments, inhibiting electronic excitation transfer to all RC pigments but 30 and 32, 
while leaving the coupling to the RP states unmodified. We found the resulting efficiency of the energy transport to 
the RP2 state in this restricted model to be indeed the same as in the original model where the excitation transport 
to the other RC pigments was allowed. This shows that only the two pigments 30 and 32 may be needed for the 
excitation energy transport to the RP states, and the role of other RC pigments appears to be mainly in the electron 
transport (described by the RP states). 


C. Amount of coherence 


According to the population current description of Equation (22), the imaginary component of coherence between 


pigments is responsible for the population currents through PSII. In the ongoing discussion about the role of coherence 
in the energy transfer in biological systems, coherence and its contribution to the energy transfer has so far been 
quantified according to different metrics that generally do not look separately at the imaginary and real parts of the 
coherence (see e.g. Refs. [SHUT] for quantification and discussion of coherence). However, as described in detail in 
Ref. [46j . the imaginary and the real part of the coherence have very different effects - the imaginary part drives 
the transfer, whereas the real part is related to the transfer-diminishing localization effects in the presence of energy 
gaps between the pigments. Therefore, we shall not only quantify the total coherence, but also look separately at its 
imaginary and real parts. 

Solving the ZOFE master equation for the full quantum dynamics gives the time-dependent electronic density 
matrix and therefore allows direct quantification of the coherences, the off-diagonal elements of the density matrix in 
the respective basis of interest. We first quantify the total amount of coherence by the sum 


c{p)= XI 


(24) 


n,m^n 


over the absolute values of all off-diagonal elements of the density matrix p, which provides an intuitive measure for 
the overall amount of coherence. This “Zi-norm of coherence” has been widely used in previous studies and is reviewed 
in Ref. [57]. 

Aside from the coherence in the site and exciton bases, we are also interested in the coherence in the domain-excition 
basis of Ref. [TB] . This is of interest since in the MRGF calculation of Ref. m, a classical rate equation is used to 
describe long-range intercomplex transfer that explicitly contains only the populations in the domain-exciton basis, 
but not the coherences. Thus, coherence in this basis is not explicitly taken into account in the calculation giving the 
population dynamics in Figure]^ (dashed lines), nevertheless good agreement with the ZOFE population dynamics 
is achieved. This is consistent with Ref. |36| . where it is discussed that the overall population dynamics, or at least 
the energy transfer efficiency, that is, the population transferred to the radical pair states, should be retained under a 
proper transition from the full quantum description to the classical rate description where coherence is not explicitly 
taken into account. 

In Figure [^, we show the total amount of coherence measured by Equaiton ( |24| ) for the site basis, the exciton 
basis, and the domain-excition basis. The domain-exciton coherence shows the coherence not (explicitly) taken into 
account in the MRGF calculations of Ref. m- 

Since the initial state is an incoherent state in the site basis, at time zero there is no coherence in the site basis. 
Similarly, in the domain-exciton basis, the coherence is initially zero. (Sites 7 and 10, which carry the initial excitation, 
belong to the same domain that consists of only these two sites.) 

In the exciton basis, on the other hand, there is initial coherence at time zero, because the excitation localized on 
sites 7 and 10 is expressed through a coherent superposition of exciton states that each have part of their excitation 
on sites 7 and 10, but are also delocalized over other sites, mainly sites of the LHCII antenna protein. 

In both site and domain-exciton basis, the amount of coherence rises initially, because the dipole-dipole interaction 
between the pigments drives the initially localized state towards a more delocalized state. During the transport dynam¬ 
ics, the total electronic coherence decreases due to interaction with the vibrations, and later due to the accumulation 
of population in the radical pair states. 
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Next, we look at the time-integrated coherences between the individual pigments and also separately at the imagi¬ 
nary and real parts of the coherence, because of their different role in the energy transfer. To this end, we time-integrate 
the absolute value of the imaginary and real parts of the elements of the density matrix in the site basis during the first 
1 ns of its propagation (at which point, as we have shown above, the majority of the excitation has been irreversibly 
trapped at RP2). In Figure]^, the resulting site-basis elements of the time-integrated density matrix are shown as 
a color matrix. The time-integrated pigment populations are the diagonal elements. One can see that there is an 
accumulation - a “damming-up” - of population on the energetically lowest sites of each protein. Such a damming-up 
effect occurs because the transport of excitation that is fast inside the complexes, where the couplings between the pig¬ 
ments are relatively strong, is slowed down between the complexes. (It also depends on the coupling to the vibrations 
and also on the gap between the “donor” protein’s pigment energies and the “acceptor” protein’s pigment energies.) 
The upper triangle in Figure]^ shows the (absolute values of) the imaginary parts of the coherences between the 
sites integrated over time. That is, it shows f |Im(p„_m/ra)| dt. As seen in Equation (22), the imaginary parts of the 
coherences give the unitary contribution to the population currents between the sites. That is, the imaginary parts 
of the coherence constitute the whole of the energy transfer, except for the transfer from the reaction center to the 
radical pair states, which is decribed by relaxation. 

The lower triangle in Figure @3 shows the real parts of the coherence, that is, / \^e[pn,m^n)\dt. We see that 
between many sites, the real part of the coherence is larger than the imaginary part. As discussed in Ref. [46], the 
real part of the coherence is related to localization effects in the energy transfer dne to the energy gaps between the 
pigments. The smaller the real part of the coherence, the smaller are these localization effects that can hinder the 
energy transfer. (Elements that are very small are ignored in the matrix in Eigurejdjl)). 

There are relatively strong imaginary parts of coherences visible in the blocks that show the coherence between 
LHCII sites and CP43 sites, and also between CP43 sites and RC sites. These are the coherences that are responsible 
for the currents between these complexes, shown by the green and red curves in Eigure j^. In Figure |4p, we see 
that there are even non-zero coherences between LHCII sites and reaction center sites. However, these are very weak 
and the resulting net population current that they cause between LHCII and the RC appears to be close to zero in 
Figure]^. 


IV. ROBUSTNESS OF ENERGY TRANSPORT AND FUNCTION IN PSII 

In the previous section we have described excitation energy transfer when initial excitation occurs in a pigment- 
protein complex on the periphery of the supercomplex and subsequent energy transfer occurs into the core complex 
via the attached LHCII trimer. In this Section, we shall explore how robust the excitation energy dynamics and the 
functional behavior of PSII are to changes in the initial state and the electron-phonon coupling. 


A. Delocalized initial excitation 

There is an ongoing discussion about the initial conditions of light absorption in light-harvesting systems, and how 
they affect the subsequent transport dynamics |36L 1391H5] . In particular, initial conditions in nature - excitation by 
sunlight - versus initial conditions in laser spectroscopy experiments - excitation through coherent laser light - have 
been discussed by a number of authors |36l ES] [48| . 

Motivated by this, we now investigate the influence of the initial state on the long-range energy transfer in PSII. 
So far, we have considered initially localized excitation of only two pigments in the LHCII, with initial coherence in 
the exciton basis but not in the site basis. Now, we contrast this by considering an initial state that is a completely 
delocalized superposition of all exciton states of all three proteins. We choose this initial state to be an (equally 
weighted) incoherent superposition of all exciton states - i.e., no initial coherence in the exciton basis. Thus, we take 
an initial density matrix, where in the exciton basis all diagonal elements (populations) are 1/A^sites (A^sites is the 
number of sites), except for the population of the radical pair states and the ground state, which are zero, and with 
all off-diagonal elements (coherences) set to zero. 

Since the exciton states are the eigenstates of the system Hamiltonian, this state would persist in the course of 
purely unitary system time propagation. But the coupling to the vibrations, the coupling to the radical pair states, 
and the radiative and non-radiative decay to the ground state cause the initial state to evolve in a non-unitary manner, 
driving the excitation transport through PSII. These dynamics lead to the creation of excitonic coherence in the course 
of time evolution, even though it is not present initially. The transport dynamics that result from this incoherent, 
delocalized initial state are summarized in Figure [6| 

The population dynamics are shown in Figure [6^ and are compared to the dynamics resulting from the localized 
initial state (Figure]^). We see that for the delocalized initial state, the population dynamics of the pigments now 
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Figure 6: As Figure]^ but the initial state is an (equally weighted) incoherent superposition of all exciton states of all three 
proteins, LHCII, CP43, and the RC. For comparison, the dashed lines in Figure A show the ZOFE dynamics from Figure [4j4. 
for the localized initial state. 


look quite different. In particular, 40% of the initial excitation is in the CP43 protein and almost 20% is in the 
reaction center pigments. As a result, 60% excitation is spatially closer to the reaction center than in the previous 
simulation. This results in a much earlier rise in the RPl population (black line) than seen in our previous simulation 
(black, dashed line). Some of this excitation, however, is initiated in higher energy states allowing for back-transfer 
from CP43 to LHCII as can be seen in Figure]^. The green curve, representing the flow of excitation from LHCII 
to CP43 is negative now for an initial period of time, showing that there is a flow of excitation from CP43 back to 
LHCII. That is to say, after initial excitation there is a net excitation transport in the opposite direction as seen 
before for the localized initial state, namely away from the reaction center. 

Figure]^ shows the time-integrated population currents between the individual pigments for the delocalized initial 
state, and should be compared to Figure for the localized initial state. The currents inside LHCII and from LHCII 
to CP43 are smaller now than before. This is expected, since now only part of the initial population is in LHCII. 
Aside from these changes, however, the overall pattern of the pathways is very similar to what we have seen for the 
localized initial state before. The currents have the same directions as before, and even their relative magnitudes seem 
very similar. This is a remarkable result, since we know in this initial condition, every pigment has an equal initial 
excitation. Nevertheless, both the short-range (intra-complex) and long-range (inter-complex) energy transfer are 
very similar for both initial conditions. Interestingly, this suggests that while there is not a single obligate pathway of 
excitation energy transport within PSH - there is a select subset of PSH pigments that are used to create an excitation 
transport network. This observation may allow future work to differentiate between pigments that are involved in 
long-range transport from those that are predominately only associated with absorbing light and then transferring 
excitation locally (and hence showing little to no net flux in either of our initial conditions). 

Despite the changes in excitation transport dynamics, the functional behavior of PSH is invariant to the change in 
initial excitation. In particular, the overall transport of population to the radical pair state RP2 is almost identical 
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Figure 7: As Figure]^ but the initial state is an (equally weighted) incoherent superposition of all exciton states of all three 
proteins, LHCII, CP43, and the RC. 


for both initial excitation conditions, as can be seen in Figure |^. (We note that because of the logarithmic time 
axis, the earlier deviation between the RPl populations for the different initial states appears to be much larger than 
the deviation between the RP2 populations at later times, but in fact these deviations are of comparable magnitude.) 
This shows that changing the initial state from the previously considered localized state to a delocalized state does 
not change the efhciency of the energy transport. The stability of the functional behavior with respect to changes in 
excitation energy dynamics can be explained by the separation of timescales between the excitation energy transport 
and the loss mechanisms. Excitation is lost in PSII on a timescale of ^ 2 ns through a combination of fluorescence 
and non-radiative decay. Excitation transport, however, occurs on a ~ 100 ps timescale. If we approximate the 
excitation transport to the radical pair states by a single rate constant, then the overall behavior is described by 
a two-rate (three-state) model - one rate for excitation to drive charge separation and a second rate of excitation 
undergoing non-radiative or radiative decay. Within this simplified model, the fraction of excitation that performs 
charge separation (as opposed to non-radiative or radiative decay) as time goes to infinity is given by 


PRP 2 {t -t oo) 


^cxc—^RP 


^exc—fRP “b ^decay 


(25) 


- the population reaching the RP2 state in the long-time limit - where fcexc-^RP is the rate of population transfer to 
the RP states and fcdecay is the rate of loss through non-radiative and radiative decay. We will use the phrase long-time 
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efficiency to refer to this quantity PRP 2 (t oo) and differentiate it from our 1 ns measure of efficiency, i.e. PRP 2 (t = 
1 ns), used in the remainder of the paper. In Figure]^ we see that at 100 ps the total excitation left on the pigments 
amounts to a population of Ppig{t = 0.1 ns) « 0.2. Inserting this value in fcexc->-RP = —yln(Ppig(f) we estimate 

the rate of excitation transport from the LHCII initial state (dashed lines) to the RP states as A:exc->-RP ~ 15 ns“^, 
which, using Eq. (25), gives a long-time efficiency of 97%. To substantially change this value, the rate of excitation 
transport to the RP states would need to slow enough to become comparable to the 0.5 ns“^ rate of excitation loss. 
For example, slowing the rate of transport to the RP states by a factor of 2 only decreases the long-time efficiency to 
93%. As a result, small modulations in the rate of transport to the RP states have only a minimal influence on the 
long-time efficiency within our model, and thus will not substantially shift the functional behavior of PSII. 

We now analyze the dependence of the amount of coherence on the initial state. In Figure [^, we see that for 
the delocalized initial state there is a significantly smaller amount of excitonic coherence involved in the transfer 
dynamics compared to the amount for the localized state in Figure (note the different range of the y axis of 
Figure]^ compared to Figure]^). This is mainly due to the fact that there is no initial excitonic coherence present 
in the case of the delocalized initial state. In Figure]^, the pattern of the time-integrated coherence between the 
individual sites appears overall similar to that for the localized initial state in Figure Ep. However, looking closer, 
there are differences, which are reflected in the large difference between the currents resulting from the delocalized and 
localized initial states. We note also, that since in Figure and D absolute values are summed up, changes of the 
signs of the individual terms that would correspond to opposite directions of the corresponding population currents 
are not visible. For example currents that lead to a flow of excitation back to LHCII in Figurej^, which have negative 
signs, cannot be identified in Figure [^, because the absolute value is taken before integrating the coherences. But 
we see in Figure]^ that just as the current between LHCII and CP43 is smaller for the delocalized initial state, so 
the imaginary parts of the coherences between LHCII and CP43 pigments are also weaker than before in Figure [4|D. 

We have also performed simulations for other initial states, including localized, delocalized superpositions with 
coherence in the site or exciton basis, and for localized initial excitations with specific phase relations that were 
previously found to lead to directed transport for models of chains of coupled sites [SS]. In all these cases, we do not 
find any significant change in the overall efficiency of the transport to the RP2 state, showing a remarkable robustness 
of the overall long-range energy transport with respect to the initial conditions. 


B. Influence of the coupling to vibrations 

We now investigate the influence of the coupling between electronic and vibrational degrees of freedom on the energy 
transfer. The coupling to the vibrations is treated approximately in our model so that it is important to consider how 
sensitive the energy transfer is with respect to changes in this part of the model. Furthermore, we expect the coupling 
to the vibrations to be an important component of the energy transfer dynamics, since it can enhance the transport, 
as is well known and has been studied for several smaller scale (single complex) systems (see e.g. Refs. [55115^ IBn| L 
In the following, we investigate these relationships by varying the strength of the coupling to the vibrations, that 
is, the magnitude of the environment spectral densities of the pigments Eq. For all pigments, we simply multiply 
their spectral densities J„(w) by a global factor that is the same for all of them, i.e., 

j;(a;)=c J„(a;), (26) 

where the constant factor c is the same for all pigments n in all proteins. As the initial state, we take again the state 
where the excitation is localized on only the two pigments 7 and 10 in LHCII that we have used in Section [lll| 

To investigate the effect of decreasing coupling to the vibrations on the energy transport, we decrease the coupling 
in stages. Figurej^shows how the energy transfer changes when we decrease the coupling to the vibrations by factors 
of ten, i.e., c = 0.1,0.01,0.001,0.0001. For each of these cases. Figure]^ shows the population of the RP2 state 
(blue, thick bar) and RPl state and the total population in each of the complexes LHCII, CP43, and the RC (colored, 
thin bars) at a time of 1 ns. We see that when we decrease the vibrational coupling, the transport efficiency - the 
population of the RP2 state at 1 ns - for a decrease of the coupling by factor 0.1 first increases slightly. Then, as 
we further decrease the coupling to the vibrations, the efficiency of the transport decreases monotonically. We can 
also see that when the vibrational coupling is smaller and smaller, there is more and more population trapped in 
LHCII. It is remarkable that even when we decrease the coupling by a factor of hundred (i.e. c = 0.01), the amount 
of population transported to RP2 within 1 ns is still higher than 60% (compared to about 80% in the case of original 
coupling strength), showing that the energy transport in PSII appears to be very robust with respect to changes 
of the coupling to the vibrations. Only when we decrease the coupling further do we find that the energy transfer 
efficiency decreases significantly. In the extreme case where we decrease the coupling to the vibrations by a factor of 
thousand (i.e. c = 0.001), it goes down to less than 20% population on RP2 after 1 ns. Decreasing the coupling further 
(c = 0.0001), causes most of the population to be trapped in LHCII and only a very small amount of population 
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Figure 8: Effect of decreasing the coupling strength to the vibrational environment by factors of ten; ZOFE simulation results. 
Figure A: Populations at 1 ns in the different complexes and of the radical pair states RPl and RP2. The bars for RP2 are 
thicker for better distinguishability. Figure B and C: As Figures [4]4, B, but the coupling to the vibrational environment is 
decreased by a factor of 10 compared to Figure [4| (corresponding to the second bar in A). For comparison, the dashed lines 
in Figure B show the ZOFE dynamics from FigurellK for the original coupling to the vibrations. Figure D and E: As B, C, 
but the coupling to the vibrational environment is decreased by a factor of 1,000 compared to the original coupling strength 
(corresponding to the fourth bar in A). In D, the dynamics for the original coupling strength are again shown as dashed lines. 


(~ 2%) reaches RP2. We can see in Figure]^ that the total length of the bars, i.e. the total population of the RP 
states and of the complexes together at 1 ns, is less than 1 (100%). The difference between the total length of the 
bars and 100% is the amount of population that is lost due to radiative and non-radiative decay to the ground state 
during the time of 1 ns. We see that this loss of population increases as the coupling to the vibrations is decreased, 
because the less population reaches the RP states, the more population remains in the excited electronic states of the 
pigments where it is subject to the decay. 














































20 



Figure 9: Transport of electronic excitation from a pigment 1 to a pigment 2 via electronic inter-pigment interaction, when 
there is a gap between the pigment’s electronic transition energies. Figure A: without coupling to vibrations, the transport 
is off-resonant - taking place between the off-resonant, purely electronic levels. This off-resonance decreases the amount of 
population transferred to pigment 2. Figure B: when there is coupling to vibrations - with adequate strength and vibrational 
energies - there can be resonance between the vibronic levels, leading to vibrationally enhanced transport. (In this sketch the 
vibrational potential surfaces are only one-dimensional, i.e., they correspond to only one vibrational mode in each pigment, 
but the same argument applies for multiple modes, i.e., more-dimensional potential surfaces. The environment spectral density 
can be decomposed into peaks of certain functional forms that represent damped molecular vibrational modes m-) 


For the two cases of c = 0.1 and c = 0.001 (the second and fourth bars from the top in Figure]^), we take a closer 
look at the energy transport dynamics in Figures]^, C (for c = 0.1) and Figures]^, E (for c = 0.001). Figure]^ 
shows how for c = 0.1 the excitation is transferred slightly faster from LHCII to CP43 and the RC than for the original 
coupling (shown as dashed lines for comparison), leading to the slightly larger amount of population reaching RP2 
within 1 ns, observed in the bar diagram Figure |8]4. The population of pigment 7, which carries half of the initial 
population, shows oscillations now, indicating that part of the population is transported back and forth between 
pigments. But this back-and-forth oscillation of a relatively small amount of population does not make the overall 
transport less efficient. Also in the current from LHCII to CP43, shown in Figure [^, there are now oscillations in 
the magnitude of the current. But the current still is always directed from the LHCII towards CP43, i.e., there is no 
backflow to LHCII, and the current is positive at all times. On the other hand, the maxima of this current are more 
than twice as large as in the case of the original coupling to the vibrations. (See Figure]^ for comparison, and note 
the different scales on the y axes). The current from CP43 to the RC and from the RC to RPl looks qualitatively 
similar to the original case. There are no strong fluctuations in the magnitude of these currents. But they are larger 
now compared to the original situation where the coupling to the vibrations is stronger. 

The transport dynamics for the case of c = 0.001 are shown in Figures [^, E. In panel D, aside from the much 
smaller amount of population transported to RP2, we can also see that it takes much longer for the population to be 
transported away from the LHCII. The reason for that becomes apparent when we look at the population of pigment 7 
in LHCII: the population oscillates back and forth between the LHCII pigments many times, and is not able to leave 
the LHCII. It is trapped in the LHCII for a quite long time due to (unitary) localization inside the complex because of 
non-uniform pigment energies. In panel E, we see that the current between LHCII and CP43 oscillates fast between 
positive and negative values. That means that excitation is flowing back and forth between these two proteins, quickly 
changing direction. The other currents are almost zero, because only a small amount of population can pass through 
CP43 and reach the reaction center. 

Why the strong decrease of the vibrational coupling leads to such a trapping, can be understood with the notion 
of vibrationally enhanced transport: the coupling to vibrations leads to transport through vibronic resonance and 
its absence leads to trapping. This concept is depicted in Figure between the excited electronic states of the 
pigments there are energy gaps that, when they are large compared to the dipole-dipole interaction strength between 
the respective pigments, cause the excitation to be trapped on one or a few pigments, because of the off-resonance of 
the electronic transitions. But if the electronic excitation of a pigment couples to vibrational modes with vibrational 
energies large enough and with a coupling strength large enough compared to the energy gaps, then the vibrational 
energy levels can fill the gaps and make the transport of excitation from one pigment to another resonant, thereby 
enhancing the transport. Furthermore, the coupling to the continuum of vibrational modes leads to relaxation down 
the energy manifold that is needed for the energy funnel, directing the transfer to the RP states, to work. 
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Figure 10: Coherence from ZOFE simulation for weaker coupling to vibrations. As Figures]^, D, but the coupling to the 
vibrational environment is decreased by a factor of 10 compared to Figure]^ 


For a delocalized initial state, in the case of very weak coupling to the vibrations, the trapping in the LHCll and 
CP43 will have smaller consequences for the tranfer efficiency, because there is also initial excitation in the reaction 
center that can be transferred to the RP states. Further simulations have shown that for a delocalized initial condition 
the population that is transferred to RP2 - in the case of the very weak coupling to the vibrations - can be three times 
higher than for the localized initial state in Figure]^. However, the population transferred to RP2 is still by more 
than a factor of two lower than when the coupling to the vibrations has the original strength. On the other hand, we 
have found in further simulations that for the original vibrational coupling c = 1, and when c = 0.1, 0.01, changing 
the initial state has little influence on the efficiency of the energy transfer. Only when the vibrational coupling is 
decreased beyond this, c = 10“^-10“"^, does the initial state have a significant impact on the overall efficiency. This 
shows again a high robustness of the energy transfer to variations of the initial condition. 

Next, let us consider how decreasing the coupling to the vibrations affects the coherence involved in the transport 
dynamics. Figure [T0| shows the magnitude of coherence for the case of c = 0.1 and should be compared to Figures]^, 
D, for the original vibrational coupling. For the weaker coupling to the vibrations (Figure [l0| ), there is more coherence 
between the sites now, because the dephasing through the vibrations is weaker. The excitonic coherence, by contrast, 
decreases (see panel A compared to Figure]^) which is expected because part of it is created by the coupling to the 
vibrations. Remarkably, these changes in the coherence do not affect the overall transport efficiency much, as we have 
seen in Figure]^ 


C. Full Markovian Lindblad simulation: larger range of coupling strength 

So far, we have only decreased the coupling to the vibrations in relation to our original model. It would also be 
instructive to investigate the energy tranfer at stronger coupling as well. Stronger coupling, however, leads to problems 
in the ZOFE calculation, because of the limited range of non-Markovianity for which the ZOFE approximation is 
valid [m [sn Eg. To still be able to gain an insight into the behavior at stronger coupling and - at least roughly 
- assess the energy transfer efficiency, we therefore use a full Markovian Lindblad calculation, where the Markov 
approximation is applied to the electron-vibration coupling. Based on the reasoning of Ref. |36j , an assessment of the 
transport efficiency can be reasonable with an appropriate Markovian description. As a check, we also calculate the 
cases of weaker vibrational coupling again with the Markovian Lindblad simulation, so that we can compare with the 
previous ZOEE simulation results. A description of the Lindblad calculation is given in Appendix]^ 

Eigure [m summarizes the effect of stronger couplings with such a Lindblad calculation, where the strength of the 
coupling to the vibrational environment is varied by a factor between 0.0001 and 50. Here, factor 1 corresponds to 
the coupling strength in the original model. We see that as before for the ZOFE calculation, the Lindblad results also 
show the highest transport efficiency in the case where the vibrational coupling is decreased by a factor of ten with 
respect to the original model. In this case, also the amount of integrated population of the RP2 state is very similar 
to the ZOFE result. However, the RP2 populations for the original coupling and for the cases of weakest coupling. 
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Figure 11: Varying the strength of the coupling to vibrational environment; fully Markovian Lindblad simulation. The strength 
of the coupling to the vibrations is varied by a factor between 0.0001 and 50. (Factor 1 corresponds to the original model.) The 
populations resulting from the Lindblad simulation are shown (in the same way as in Figure]^). The initial state is the same 
localized initial state as for Figure where all population is on the two pigments 7 and 10 in the LHCII, without coherence in 
the site basis. 


calculated with Lindblad, are in less good agreement with the ZOFE results. 

Let us now look at the cases where the vibrational coupling is increased compared to the original model. In 
Figure [TH we see that as the coupling is increased, the transport efhciency becomes lower and lower, until for a fifty¬ 
fold increase, it is almost zero. At the same time the amount of population that stays in the LHCII becomes larger and 
larger. This suppression of the transport when the coupling to the vibrations is stronger may be understood in terms 
of the quantum Zeno effect |53j . where the electronic excitation interacts with the vibrations - i.e., it is “measured” 
by the vibrations in the site basis - so strongly, that is, at such a high rate, that the excitation is forced to stay in 
the original state and cannot move anymore. In terms of the expression for the population currents. Equation (221, 
this effect can be understood from the fact that the stronger coupling to the vibrations causes stronger dephasing, 
destroying the coherence between the sites that drives the energy transfer, and thereby diminishes the energy transfer. 
Because of this effect at stronger coupling and the trapping of the excitation at weak coupling, the highest transport 
efficiency is obtained in the intermediate regime. There, the magnitude of the coupling to the vibrations is such that it 
gives the optimal balance between vibrationally enhanced transport (as described in Figure]^ and the quantum Zeno 
effect. These results of the Lindblad calculation therefore provide more evidence that the strength of the coupling to 
the vibrations that we used in our original model lie roughly in the region of optimal coupling strength. However, 
remarkably, these results again show that the energy transfer efficiency could be even (slightly) higher at a factor 
~ 10 weaker coupling to the vibrations. 


V. SUMMARY AND CONCLUSIONS 

In this work, we simulated the long-range inter-complex energy transfer in the PSH supercomplex - from LHCII 
via the core complex CP43 to the reaction center - by means of a full quantum calculation, using a non-Markovian 
ZOFE quantum master equation description. We found that nearly all energy from the initial electronic excitation is 
transported to the reaction center within about 200 ps, and that after 1 ns, about 80% of the energy is transformed 
into charge separation in the reaction center. Less than 5% is lost through radiative and non-radiative decay within 
this time. 

These findings for the overall energy transfer dynamics are in good agreement with the results of Ref. m, where a 
modified-Redfield-generalized-Forster (MRGF) rate-equation method was applied to calculate the long-range energy 
transfer in the entire PSH supercomplex. In the MRGF method, strongly coupled pigments are assigned to domains 
treated with a modified-Redfield description, whereas excitation transfer between the domains, where the couplings are 
weaker, is treated with a generalized-Forster description. In the present work, to compare our ZOFE quantum-master- 
equation results to the MRGF rate-equation description of Ref. [TS] , we used the MRGF method to calculate the energy 
transfer in the truncated supercomplex LHCH-CP43-RC of PSH that we considered. The resulting energy transfer 
dynamics are in very good quantitative agreement with that calculated with the ZOFE master equation, showing 
that the simpler MRGF rate equation provides an adequate description of the quantum energy transfer dynamics 
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and energy transfer efficiency of PSII. The observation that the MRGF results agree well with the ZOFE results 
shows that the subdivision into domains, which is related to the notion of supertransfer, i.e., coherent delocalization 
inside domains leads to enhanced transfer rates between the domains [3S] , is a reasonable approach to describe energy 
transfer in large light-harvesting complexes, where the coupling between the pigments can be very inhomogeneous. 
We have further shown that excitation energy transfer is significantly slower within a pure Forster calculation, which 
further supports the importance of supertransfer for rapid excitation transport within PSIL The close agreement 
between ZOFE and MRGF is also consistent with the hndings of Ref. [35], where a full quantum description of 
electronic excitation transfer is compared to that of a classical rate equation, and it is shown that the rate equation 
description is capable of capturing the overall features of the energy transfer when properly derived from the full 
quantum description. 

In addition to the commonly used description in terms of time-dependent populations of the excited electronic states 
of the pigments, we also analyzed the energy transfer dynamics in terms of excitation population currents between the 
pigments. This analysis revealed the pathways and directions of the energy transfer between the individual pigments 
and the net fow of energy between the different complexes. We found that when initially all excitation is in the 
LHGII antenna complex, the net flow of energy is directed towards the reaction center at all times and no temporary 
intercomplex backflow occurs. Integrating the population currents over time showed that there is a complex network 
of pathways rather than just one or two dominant pathways, in accordance with the findings of Ref. |16j . The energy 
transfer between different proteins was found to rely on several different pathways. Interestingly, instead of being 
uniformly directed towards the reaction center, often excitation energy is transported in loops between pigments, i.e., 
excitation is transferred from one pigment to a neighboring pigment and from there to the next, etc., in such a way 
that it goes around in a closed circuit of pigments and in the end comes back to the initial pigment. Such loops were 
found in all three subcomplexes that we considered, LHGII, GP43, and the reaction center. In the reaction center, 
the efficient excitation energy transport to the charge separated states appeared to rely mainly on only two pigments 
(30 and 32), providing evidence that the rest of the reaction center pigments are needed primarily for their role in the 
electron transport rather than for the excitation transfer; some of these reaction center pigments appear to constitute a 
strong loop excitation current, but do not contribute to the net transport of excitation energy to the charge separated 
states. To test this conclusion, which was drawn from the analysis of the pathways of integrated population currents, 
we ran additional simulations with all reaction center pigments - except for 30 and 32 - decoupled from the rest 
of the pigments so that transfer of excitation to and from these decoupled pigments was inhibited, but leaving the 
coupling to the charge separated states unmodihed. The resulting efficiency of the energy transport to the charge 
separated states in this restricted model was found to be the same as in the original model where excitation transport 
to the other reaction center pigments was allowed, showing that the excitation energy transport to and inside the 
reaction center appears to rely indeed only on the two pigments 30 and 32. This confirmed the conclusions drawn 
from the analysis of the integrated population currents and showed that this analysis is suited to identify transport 
pathways and assess their importance, providing a useful tool for the investigation of design-function relationships in 
light-harvesting systems. 

By simulating the energy transfer with the ZOFE master equation, we were able to quantify the amount of electronic 
coherence that is involved in the energy transfer dynamics. Motivated by the description of the energy transfer in terms 
of population currents and the different roles of imaginary and real part of the coherence presented in a companion 
paper [46j . we quantified the real and imaginary components of the coherence between the pigments separately. 
We found that the real components, which are related to localization of excitation due to energy gaps between the 
pigments |46j . are overall larger than the imaginary components, which are proportional to the population currents 
between the pigments and drive the excitation energy transfer in the PSII supercomplex. 

To find out how the energy transfer in PSII depends on the initial condition - which is determined by the specihc 
features of the light used for the excitation - we ran simulations of the energy transfer for different initial states; first 
we considered initial excitation that is localized on only two LHGII pigments; then we started from an initial state in 
which the excitation is completely delocalized over the LHGII, CP43, and reaction center pigments. We found that 
even though for the delocalized initial state the energy transfer dynamics of course look different from that for the 
localized initial state, the overall efficiency of the energy transfer to the charge separation in the reaction center is 
nevertheless the same in both cases. Also the overall pattern of the energy transfer pathways between the individual 
pigments is very similar for the different initial conditions ~ for both short-range (intra-complex) and long-range 
(inter-complex) energy transfer. This shows that the energy transfer in PSII is remarkably robust with respect to the 
initial conditions. 

To learn how sensitive the energy transfer dynamics are to variations in the coupling between the electronic and 
vibrational degrees of freedom, we varied the strength of this coupling to the vibrations in our simulations: decreasing 
the coupling to the vibrations for all pigments by a factor of ten compared to the original situation did not change 
the energy transfer efficiency much - it even slightly increased the efficiency, showing that the energy transfer is very 
robust to variations of the coupling to the vibrations. Even when the coupling to the vibrations was decreased by a 
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factor of 100, the amount of energy transferred to the charge separation within 1 ns decreased only by about 10%. 
We assign this to a decrease of the well known effect of vibrationally enhanced transport, which helps the transfer 
of energy through more resonant energy levels when coupling to vibrations (vibrational energy levels) is present. 
Only when the electron-vibration coupling was decreased by a factor of 1,000, does the amount of energy transferred 
decrease significantly (by up to a factor 4 compared to the original situation), due to trapping of the excitation in 
the LHCII complex. This trapping is characteristic of the predominantly unitary dynamics occuring in the absence of 
vibrationally enhanced transport. These observations show that for PSII the strength of the coupling to the vibrations 
that is needed for the effect of vibrational enhanced transport to direct the energy transfer efficiently to the reaction 
center is rather small. 

However, in contrast to the Fenna-Matthews-Olson (FMO) complex, where calculations have revealed that the 
transport efficiency does not fall under 70% when the coupling to the vibrational environment is decreased, even to 
extremely small coupling strength [40] , for PSII we found that the efficiency can become very low - almost zero - for 
extremely small coupling to the vibrations. This shows that contrary to FMO, for PSII the coupling to the vibrations 
is crucial for achieving efficient transport. 

To assess how the energy transfer changes when the coupling to the vibrations is increased compared to the original 
situation, we performed Markovian Lindblad simulations, instead of the ZOFE calculations because the latter were 
not applicable in this parameter range. Here, we found that if the coupling to the vibrations is stronger than in the 
original situation and is further increased step by step, the efficiency of the energy transfer decreases monotonically, 
until for a fifty-fold increase of the coupling, almost no energy is transferred to the reaction center anymore (within 
the reference time of 1 ns) due to the strong dephasing that destroys the coherence needed for the energy transfer [45] . 
These approximate Lindblad simulations, together with the ZOFE simulations, thus provide evidence that the original 
strength of the coupling to the vibrations lies in the parameter region for optimal energy transfer efficiency. We showed, 
however, that for a coupling that is about ten times weaker than in the original situation, the efficiency could be even 
slightly higher. This could be an important feature to reflect upon, in the light of the discussion about the possible 
(evolutionary) optimization of natural light-harvesting systems like PSII. 


Appendix A: Pure Lindblad simulation of excitation energy transfer 


To be able to roughly assess the efficiency of excitation energy transfer in our truncated PSII model for an extended 
range of the electron-vibration coupling, in Section jlV C| we used a pure Lindblad simulation for the excitation 
energy transfer, applying the Markov approximation to the electron-vibration coupling. Here, we briefly describe the 
calculation and how the corresponding parameters are obtained. 


For the simulation, we employ a Lindblad master equation (19), where we now include Lindblad terms describing 


the coupling of the electronic excitation to the vibrations of the pigments and protein environment - in addition to 
the Lindblad terms for transfer to the RPl and RP2 states and radiative and non-radiative decay to the ground state 
(relaxation terms) that were also included in the ZOFE simulations. The new Lindblad terms describing the electron- 
vibrational coupling contain Lindblad operators = \ n){n \ and coupling parameters for each pigment n, so 

that - instead of Eq. (20) for the ZOFE simulation - we now have the pure Lindblad master equation 


/ 1 1 \ 
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where the terms in the second line again describe transfer to the RPl and RP2 states and radiative and non-radiative 
decay. In the Lindblad term in the first line, the electronic excitation of each pigment n is coupled to vibrations with 
respective coupling strength To calculate the excitation energy transfer from the full Lindblad description (Al), 

we need to assume values for the coupling parameters In the following, we show how we obtain these parameters 

through a simple approximation from the parameters that we used to describe the non-Markovian electron-vibration 
coupling within the ZOFE description. To this end, let us again consider the special form of the environment 
correlation function (16) used in the ZOFE description 
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The pure Lindblad description (Al) is obtained in the Markov limit, where 

an{t-s) = 7)(“‘^(5(t- s), 


(A3) 
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i.e. an{t — s) decays fas t com pared to the timescales of the system dynamics. To express the parameters 7 ^“'^ through 
the parameters of Eq. (A2| in this limit, we consider the Dirac series Se{x) = ^ exp(—|a:|/e) for e —)• 0, which lets us 
write 




Inserting Eq. (A4) in Eq. (A2) we have 
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Inserting (A5| in the Eq. (14) for the ZOFE auxiliary operator, we get 
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since = Ln (see Eq. (15)). Since we know that in the Markov limit a„(t — s) = — s), the auxiliary 
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operator of the ZOEE mas ter e quation becomes the constant operator Oq”^ (t) = (see Eqs. (14) and (15)) [47 


from comparison with Eq. (A6) it follows that we can approximate the needed coupling parameters forT 
descriptions as 
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i.e. in the case where the 7 „j are large compared to the system energies. We note that in order for the 7 „ to stay 
finite in the 7 „j —>■ oo limit, not only the 7 „j, but also at least one of the parameters has to go to infinity. Inserting 


the parameters for the electron-vibration coupling of the ZOEE simulation (according to Section IIB2) into Ec 
we find the following values for the parameters that we used for our full Lindblad simulation in Section 
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4, 200 cm“^ for the Chi A pigments in LHCII, 

5, 280 cm“^ for the Chi B pigments in LHCII, 

7 ""‘"“ « 140 cm“^ for the pigments in CP43, 

and « 190 cm“^ for the pigments in the reaction center. 

We can see from the fact that these values are not at all large compared to the system energies (the differences of 
the eigenenergies of the system Hamiltonian are in the range of up to 1,300 cm“^) that the Markov approximation 
is not well justified here and that therefore our full Lindblad simulations should be considered to be a rather crude 
approximation of the non-Markovian dynamics. 


Appendix B: Effect of coupling to the radical pair states 

To see the influence of the radical pair states on the transport dynamics, we now look at the dynamics with the 
coupling to the RP states switched off. That is, there is no transfer of population to the radical pair states. 

Figure [T^ shows the resulting dynamics for the same quantities as in Figure [^ but with the coupling to the RP 
states switched off. We see in Figure [T2j4 that up to a certain time, the populations in the proteins with and without 
coupling are equal. The further a protein is away from the RP states, the later the population without the coupling 
deviates from the population where the coupling is present. Instead of falling off fast, now - without the coupling to 
the RP states - after around 200 ps the populations decline much slower; this residual, slower decline is due to the 
radiative and non-radiative decay of the pigment excited states. After 200 ps, the residual, average population per 
pigment is highest in the RC, and it is lowest in LHCII. This is because the average energy in the RC is lowest and 
that in LHCII is highest, providing the “energy funnel” that directs the transport to the RC. 

Let us next look at the coherence. When we compare the total coherence in Figure |T2p, where the coupling to the 
RP states is switched off, with that in Figure]^, where the coupling to the RP states is present, we see that up to 
about 10 ps the coherence with and without coupling is equal. After that, however, the coherence without coupling 
to the RP states goes to finite values that decrease only slowly, whereas in the case when the coupling is present, the 
coherence falls off to zero. This behavior can be explained with the Cauchy-Schwarz inequality 


(Bl) 

















26 


1.0 


.2 

’■+J 

3 

Oh 


0.0 


A" 


..... . . 


- LHCIIC 



- CP43 

\ y 


- RC 



- RPl 

\ ^ \ - 


- RP2 

- LHCII site 7 

\ / 
\ ' 




.1 


■-1^- 


0 10 - 


10 ° 


10 ^ 
time [ps] 


10 ^ 


10'^ 


O 

o; 

a; 

o 

o 

3 

o 



time [ps] 


Figure 12: As Figure and C, but no coupling to the RPl and RP2 states. For comparison, the dashed lines in Figure A 
show the ZOFE dynamics from Figure |4|\ (solid lines in Figure |4]4) where the coupling to the RP states is present. 


that holds for the coherences Pij{t) and the populations puif), Pjj(t) of any density matrix p{t). It says that coherences 
can only ever arise between states that are populated. The larger the populations of two given states are, the larger 
can the coherence between these states potentially be. In the case where there is transfer to the RP states, population 
is taken away from the excited states of the pigments by the RP states, and consequently, the coherence, being limited 
by the population through the Cauchy-Schwarz inequality, has to fall off to zero. On the other hand, when the transfer 
to the RP states is switched off, the population remains in the excited states of the pigments and the coherence stays 
finite. 
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